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SELF-CALIBRATING METHOD FOR PHASE ESTIMATION FROM 
MULTIPLE PHASE-SHIFTED PATTERNS 

Technical Field of the Invention 

The present invention relates to phase demodulation, and in particular to phase 
demodulation of a sequence of fringe patterns with arbitrary and unknown phase shifts 
between them. Such patterns typically occur in optical interferograms, but can also occur 
5 in interferograms produced by more general waves such as acoustic, synthetic aperture 
radar (SAR), ultrasound, and x-ray. More general image sequences, such as video, can in 
some cases also be modelled as fringe pattern sequences. 

Background Art 

10 Interferometry refers to experiments involving the interference of two light 

waves that have suffered refraction, diffraction or reflection in the interferometer. Such 
experiments typically involve the testing of lenses, mirrors and other optical components. 
The principle parameter of interest in interferometry is the phase difference between 
interfering beams. 

15 Present methods of demodulating sequences of phase-related fringe-patterns 

generally rely on a priori knowledge of the phase-shifts between each of the individual 
patterns. When these relative phase-shifts are available, it is possible to estimate the 
spatial phase by using a generalised phase-shifting algorithm (PSA). In cases where the 
relative phase-shifts are not a priori known, the estimated spatial phase will be incorrect 

20 unless special error-reducing PSAs are used or methods to estimate the actual phase-shift 
between patterns are used. 

Error-reducing PSAs are useful where the deviation of the actual phase-shift 
from the expected phase-shift is small (typically less than 0. 1 radian) and deterministic. 
When the phase-shift deviation is larger and/or non-deterministic (ie essentially random), 

25 then other methods must be used. 

Methods have been developed to estimate phase shifts between phase-related 
fringe patterns. Certain of the methods, based on statistical (least squares) estimation, 
requires at least 5 fringe patterns to work. Methods based on fitting the phase shifts to an 
ellipse also require at least 5 fringe patterns. Methods based on the statistical properties 
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of spatial phase histograms have been proposed, but are heavily dependent on certain 
restrictive, and often unrealistic criteria. More recently methods based upon the 
correlation between patterns have been proposed but are very sensitive to the number of 
full fringes in the frames. These correlation methods assume that the fringes are 
5 essentially linear in the frame, which means only trivial fringe patterns can be processed. 
Methods based on the maximum and minimum variation of individual pixels in different 
patterns typically require at least 1 5 frames to operate effectively. 

In summary, current methods are unsatisfactory in many situations, especially for 
a small number of phase-related fringe-patterns, and many algorithms do not work in the 
1 0 case of three or four frames. 

Disclosure of the Invention 

It is an object of the present invention to substantially overcome, or at least 
ameliorate, one or more disadvantages of existing arrangements. 
15 According to a first aspect of the invention, there is provided a method of 

estimating relative phase shifts between fringe pattern images in a sequence of phase- 
related fringe patterns, said method comprising the steps of: 

a) removing offsets from each of said fringe pattern images to obtain pure 
AMFM patterns; 

20 b) determining contingent analytic images corresponding to each of said 

AMFM patterns; 

c) determining phase differences from dependent pairs of said contingent 
analytic images; and 

d) estimating phase shifts between pairs of said contingent analytic images, 
25 wherein said phase shifts between pairs of said contingent analytic images are said 

relative phase shifts between fringe pattern images. 

According to a second aspect of the invention, there is provided a method of 
estimating a spatial phase of fringe pattern images in a sequence of fringe patterns, said 
method comprising the steps of: 
30 a) removing offsets from each of said fringe pattern images to obtain pure 

AMFM patterns; 
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b) determining contingent analytic images corresponding to each of said 
AMFM patterns; . 

c) determining phase differences from dependent pairs of said contingent 
analytic images; 

5 d) estimating phase shifts between pairs of said contingent analytic images; 

and 

e) estimating a spatial phase of said fringe pattern images. 
According to another aspect of the invention, there is provided an apparatus for 
implementing any one of the aforementioned methods. 
10 According to another aspect of the invention there is provided a computer 

program product including a computer readable medium having recorded thereon a 
computer program for implementing any one of the methods described above. 

Brief Description of the Drawings 

15 A preferred embodiment of the present invention will now be described with 

reference to the drawings, in which: 

Fig. 1 shows the local structure of a typical fringe pattern; 
Fig. 2 shows the Fourier domain coordinates; 

Figs. 3A and 3B show the phase shift parameters 8 n for the case where N=5 on a 
20 polar plane; 

Fig. 4 shows a flowchart of an automatic calibration method in accordance with a 
preferred embodiment; 

Fig. 5 shows an example of a fringe pattern; 

Fig. 6A and 6B show the magnitude and phase components respectively of an 
25 analytic image of the example fringe pattern shown in Fig. 5; 

Fig. 7A shows a phase difference between dependent analytic images; 

Fig. 7B shows a histogram of a spread in values of the phase difference shown in 

Fig 7A; 

Fig. 8A shows a modulus of the phase difference between dependent analytic 
30 images shown in Fig. 7A; 

Fig. 8B shows a histogram of a spread in values of the modulus of the phase 
difference shown in Fig. 8B; 
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Fig. 9 shows a polarity function calculated from the phase difference shown in 
Fig. 7A and the modulus of the phase difference shown in Fig. 8 A; 

Figs. 10A and 10B show estimations of an amplitude modulation and a phase of 
the fringe pattern shown in Fig. 5; 
5 Fig. 1 1 shows an estimate of an offset of the fringe pattern shown in Fig. 5; 

Fig. 12 shows a simple weight function calculated from the estimated errors in 
the analytic image shown in Figs. 6 A and 6B; and 

Fig. 13 shows a schematic block diagram of a general purpose computer upon 
which the preferred embodiment of the present invention can be practiced. 

v 10 

Detailed Description including Best Mode 

Interferometry is the use of interference phenomena for measurement purposes. 
This allows for the measurement of very small angles or small displacements of objects 
relative to one another. An interferometer is a device to make such measurements, which 

15 typically uses a beam of light, coming from a single source, such as a star, a laser or a 
lamp, and two or more flat mirrors to split off different light beams. These beams are 
then combined so as to 'interfere' with each other. The beams interfere with each other 
by constructively adding together and cancelling each other out, thereby forming 
alternating bands of light and dark, called fringe patterns. 

20 What makes the interferometer such a precise measuring instrument is that these 

fringes are only one light-wavelength apart, which for visible light corresponds to about 
590 nm. With the fringe patterns having a cyclic nature with phases that shifted, the 
important measurement parameter is the relative phase-shifts between these phase-related 
fringe patterns. 

25 Consider an intensity f n (x 9 y) of the n th fringe pattern in a sequence of phase- 

related fringe patterns with the following mathematical form: 

f n (x,y) = a{x 9 y) + b(x 9 y)cos[y/(x 9 y) + S n ] (1) 

wherein a(x 9 y) is a background or offset, b(x,y) is the amplitude modulation, and 
y/{x 9 y) is the phase. Each fringe pattern has a phase shift parameter S n . An example of 
30 such a fringe pattern is shown in Fig. 5. 
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A improved method of estimating relative phase-shifts between phase-related 
fringe patterns f n is proposed. The method is based on an isotropic quadrature transform. 
The isotropic quadrature transform is used to estimate the spatial phase [y/(x,y) + 8 n ] in 
one fringe pattern f n , followed by comparing the spatial phase [y/(x,y)+S n ] between 
fringe patterns f n to estimate the relative phase-shift(s). With 3 or more frames, the 
method is independent of background variations a(x 9 y) . 

The exact value of phase y/(x,y) can be calculated from three or more fringe 
patterns/! if the values of the phase shift parameter S n is known for each fringe pattern f n . 
For more than three fringe patterns f n (x,y) , the phase y/(x y y) is over-determined and is 
usually estimated in a maximum likelihood process. Other estimation processes may be 
more appropriate in cases where additional systematic error or distortion effects apply to 
Equation (1). 

Consider a least squares estimate for the phase y/(x,y) from N fringe patterns: 
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where c n = cosJ w , and s n =sinJ„, and the summations are over N fringe patterns ie. 
^ = - Least squares Equation (2) can be inverted, providing c n and s n are such 

that a singular matrix does not occur, as follows: 
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Equation (3) allows for the offset a(x 9 y), the amplitude modulation b(x y y), 
and the phase y/(x y y) to be calculated. However, before Equation (3) can be used, all the 
phase shift parameters S n must be known. 

The phase shift parameter 5 n of each fringe pattern is estimated by using a 2-D 
quadrature technique based on a spiral phase Fourier operator. Estimation of phase 
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iff(x,y) in the above manner has many advantages related to the extraordinary 
demodulation properties of the spiral phase Fourier operator. 

A first step in determining the phase shift parameters S n is to remove the 

background (offset) a (x, y) from each of the fringe patterns f„(x,y). In the preferred 
5 embodiment, this is done by forming inter- frame differences between pairs of fringe 
patterns f n (x y y) as follows: 

g nm =/,-/*= ~g mn = b {cos [y/ + S n ] - cos [y/ + S m }} (4) 

The inter-frame difference function g nm is a pure AMFM (amplitude modulation 
and frequency modulation) function without offset and therefore can be demodulated if 
10 the inter-frame difference function g nm , which is also a fringe pattern, is locally simple. 

For the fringe pattern g nm to be locally simple, it has to have a unique orientation fi y) 
and spacing X{x,y) at each location (x,y), except, perhaps, at a finite number of 
locations where the orientation fi(x,y) and/or spacing Pi{x y y) may be uncertain. 
Advantageously, local simplicity is a property which distinguishes fringe patterns from 
15 other, more general, patterns. Fig. 1 shows the local structure of a typical fringe pattern 
indicating the orientation (3{x,y) and spacing Z(x,y) for the small region. 

A spiral phase Fourier operator V{}, hereinafter referred to as a vortex operator, 
is defined as follows: 

V{f(x,y)} = ¥- 1 {exp[i<p]F{f{x,y)}} 4 
F{f{x>y)} = F ( u > v )= J j f(x,y)exp[-27Vi(ux + vy)]dxdy 

F" 1 {F(u 9 v)} = f(x 9 y)= J J F(u 9 v)expl+2xi(ux + yy)]dudv 

u — qcos<j) 
v — qs\r\(/) 

20 The Fourier operator F { } transforms a general function f(x,y) from the spatial 

domain (x,y) to the spatial frequency domain (m,v). The Inverse Fourier operator 
F" 1 { } transforms a spatial frequency domain signal back to the spatial domain (x,y) . In 




504654.doc 



-7- 



practice, for discretely sampled images, the Fourier and the Inverse Fourier transforms are 
efficiently performed by a Fast Fourier transform algorithm and an Inverse Fast Fourier 
transform algorithm. 

The spiral phase exp[/0] is wholly a function of the angular component of the 

5 spatial frequency polar coordinates (q,<p) shown in Fig. 2. The vortex operator V{} is 
applied to the inter-frame difference functions g nm to obtain: 

V {g mn } = i exp [ifl] b [sin ty + 8 n )- sin (y/ + 5 m )] (6) 

From Equation (6) it can be seen that the vortex operator V{} is an isotropic 
quadrature operator in that it converts cosine fringes at any orientation J3(x,y) into sine 

10 fringes at the same orientation J3(x 9 y). The approximation in Equation (6) is valid 
except in regions where the amplitude modulation b(x,y) , spacing A(x,y) or orientation 
P ( x >y) var Y to ° rapidly. Areas where the radius of curvature of the fringe pattern g nm is 
smaller than the fringe spacing X{x,y) also give a lower accuracy approximation. 
However, for most typical fringe patterns g nm , the approximation is good and any errors 

15 are isotropically distributed. In contrast, conventional methods, such as the half-plane 

Hilbert transform, cause highly directional errors, and are not suitable for high precision 

demodulation of general fringe patterns. 

Next, the local orientation /? {x*y) is compensating for. To do so, the orientation 

/3(x,y) must be estimated. There are a number of ways to do this. One method relies on 
20 the gradient of the fringe pattern g nm . The ratio of the x and y components gives the 

tangent angle fi{x,y) from the local Taylor's series expansion of the phase y/(x,y) 

about the point (x 0 , y 0 ) : 

y/(x,y) = y/ 00 +2?vu 0 (x-x 0 ) + 2ttv 0 (y-y 0 ) + O(r 2 ) 

dV '^ ,2nu 0 (7) 



dx 
dy/(x,y) 



y=y 0 



dy 



= 2tcv 0 



X=X 0 
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The gradient components are: 

^ S b(-™(v + 6 a ) + 8m(y + 6 m ))^ 

The ratio of the gradient components gives the tangent angle: 

%^/%^ = ^ = tan/? nm (9) 
dy I dx u 0 

Again the approximation of Equation (9) is valid in regions with smoothly 
varying parameters as defined for the vortex operator V{}. Other methods of orientation 
estimation may be used. For example, methods based on the multi-dimensional energy 
operator or the vortex operator itself may be used. Generic orientation estimation is 
assumed at this point, and the orientation /? for any fringe pattern f n or inter-frame 

difference function g nm may be calculated. Alternatively the orientation /? may be 
estimated from combinations of the previous methods in a statistically advantage manner. 
It will be assumed hereinafter that there is just one orientation estimate, which shall be 
called orientation estimate f3 e . The orientation estimate fi e typically contains an 
ambiguity in the sign: 

exp [iJ3] = cos J3 + / sin J3 = ± U ° +IV » ( 1 0) 

V"o +v o 

Different orientation estimates f3 e typically have the sign ambigufties in different 
places. It is worth noting that the overall sign ambiguity arises from the general difficulty 
in distinguishing a fringe at orientation (3 from a fringe at orientation fi + iz . 

Equation (4) is rewritten as: 

g nm =26sin[> + +S„)/2]sin[{S m -S n )/2] (1 1} 

and Equation (6) as: 

exp[-^JV{g„J = 2/6exp[/(^->9j]cos[^ + (^ + J n )/2]sin[(^-J (1 )/2] (12) 
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The complex exponential in Equation (12) has possible values depending 
on the orientation estimate (3 e , and shall be called a polarity factor h(x,y) : 

h{x,y) = oxp[i{j3-p e )] (13) 

A contingent analytic image g nm is defined as: 

5 g nm =i(g nm -exp[-ifl e }/{g nm }) (14) 

By substituting Equations (11), (12) and (13) into Equation (14), the contingent 
analytic image g nm is rewritten as: 

g nm = 2bsm[{S m -5 n )l2\hcos[y, + {S m + S n )/2]+isin[ V , + (S m +S n )/2]} (15) 

This definition of an analytic image is contingent upon the orientation function 
10 J3 e used in the definition shown in Equation (14). It is noted that the contingent analytic 
image g nm can be achieved by placing the inter-frame difference function g nm in the 
imaginary part of a complex image, and placing the imaginary output of the vortex 
operator V{} in the real part. 

A phase a nm of the contingent analytic image g nm is calculated as follows: 

15 a nm =Arg{ig nm -i e xp[-ij3 e ]v{g nm }} = h[^ + {S m +S n )/2'] (16) 

The calculation of phase a nm is possible as long as sin[(£ m -S a )/2^ ^ 0 , which 
only occurs if S m -5 n * 2k j where j is an integer. Clearly the failure to calculate phase 

4 

a nm in such a case is due to the phase shifts S m and 8 m being essentially equal and the 
inter- frame difference function g nm being identically zero. 
20 The phase difference between dependent pairs of contingent analytic images, 

such as g nm and g mk , is calculated as follows: 

cc nm -a mk =h[ W + {S m + 5„ )/2] -h[y+(5 t + 8 m )/2] = h [(S n - 5 k )/2] ( 1 7) 

This allows a difference in phase between any two fringe patterns f n and fk to be 
found, namely (8 n -5k), with a sign ambiguity as a result of the polarity factor h. It is noted 
25 that all phase calculations, including phase differences, are calculated modulo 7.K . 
Various methods may be used to extract both the polarity factor h and the phase 
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differences (S n -8k). In a simplest embodiment, from Equation (17), the absolute values of 
the phase difference between two fringe patterns f n and fa are averaged over the full fringe 
pattern area S: 

\5 n -S k \ =2 JJsl rr (18) 

I k\mean JJ ^ 

5 A preferred embodiment uses a weighted integral with weighting w(x,y) over the 

fringe pattern area S. The weighting w(x,y) used should be a measure of the estimated 
accuracy of the local estimate of the phase a nm . This is related to the AMFM 
characteristics of the fringe patterns f„(x,y) in the sequence of phase-related fringe 
patterns. Thus, the absolute values of the phase difference {S n -S k ) between two fringe 
10 patterns/, and f k > weighted averaged over the full fringe pattern area S is: 

\°n ~ °k Righted ~ I rf ~ — ~ U * ) 

mean })s w xX'yjdxdy 

This allows for a very low (or zero) weighting w(x,y) to be used where the phase 
a nm varies quickly, as is the case for discontinuities, or areas where the modulation 
b{x,y) is low. Equation (17) can be rewritten as: 

15 a nm -a mk , ( ^Az^ (20) 

By arbitrarily taking the phase difference (5 n -Sk) for these particular values of // 
and k as positive, from Equation (20) follows an estimated polarity factor h e as: 

h e ^y) = ^ m ~ am \ (21) 



This estimated polarity factor function h e can be used in all further computations, 
20 because a global sign convention has effectively been set for the phase difference (S n -S k ) 
in the fringe pattern sequence. 

Next, the phase difference 8 n - S k may be calculated for all possible pairs. The 

minimum requirement is to calculate the phase difference S n -S k for all sequential pairs, 
which amounts to N values for N fringe patterns f n . By also calculating the phase 
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difference 8 n —S k for non-sequential pairs, additional statistical robustness in the 
estimates is achieved. The number of combinations of two-phase differences 8 n - 8 k 
from N phases are: 



In the particularly interesting case where N=3, the number of combination is also 
3, which equals the number of sequential differences too. Figs. 3A and 3B show the 
phase shift parameters S n for the case where N=5 on a polar plane. The points 1-5 are 

thus on a unity circle at points corresponding. to exp(iS n ). The connecting line between 
points n and m corresponds to [zxp(iS m )-exp(iS n )]. The length of each connecting line is 
proportional to sin[(£ m - S n )/l~j, which also appeared in the inter- frame difference 
function g nm of Equation (1 1). The length of each connecting line is a direct indication of 
the reliability of the estimate of the phase difference 8 n —8 k , as the length is the factor 
appearing in the magnitude of the inter-frame difference function g nm . Fig. 3 A shows the 
connecting lines for all sequential pairs, whereas Fig. 3B shows all possible connecting 
lines, and thus inter- frame differences 8 n -8 k , which adds another 5 phase differences 

8 n — S k . With the phase differences y mn defined as y mn = S tn -S lt , an estimation y mn of 

each of the phase differences y mn may be calculated using a weighted mean, *and 
compensated by the estimated polarity h e as follows: 



In situations where the fringe patterns f n contain noise and/or other distortions, 
the best phase differences estimates y mn may come from a two dimensional fit of all the 
points corresponding to exp(/£„), as shown in Figs. 3A and 3B, to a unit radius circle. 

Without loss of generality, one of the phase shift parameters 8 n may be set to 
zero. For convenience the first phase shift parameter 8 X = 0 is set. This is consistent with 
the idea that the absolute value of the phase y/(x,y) is undefined and only phase 



2(N-2)\ 



(22) 




4 



(24) 
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differences y mn —8 m —8 n are really meaningful. The remaining (relative) phase shift 
parameters 8 n are calculated as follows: 

<?2 = ^1 + f 21 

S z ^S 2 +f 22 (25) 
5 A = S 3 + f 43 

These values of the phase shift parameters <^are substituted back into the 
5 general PSA defined by Equation (3). This allows for all the fringe pattern parameters, 
namely the offset a(x,y) , the amplitude modulation b(x,y) , and a phase y/'(x,y) to be 
calculated. It is noted that the calculated phase y/\x,y) then includes a common phase 

shift, namely the phase shift Sj. 

However, there may be some residual calibration error due to the original vortex 
10 operator failing in certain areas. These errors are rectified by re-evaluating the weight 
function w(x,y) based on the present estimates of the offset a(x,y), the amplitude 

modulation b(x,y) , and phase y/'(x,y) , as well as their partial derivatives. This allows 
for the assessment of how well the original fringe patterns, in the form of the inter-frame 
difference functions g nm , fulfilled the smoothness requirements of the vortex operator 

15 V{}. Areas where the phase a nm varies too quickly, or areas where the modulation 
b(x,y) is too low, the weighting function w(x,y) is reduced or set to ze#o. Once a new 
weight function w(x,y) is formed, the phase calibration process, from Equation (19) 
onwards, is repeated, resulting in better estimates for the phase shift parameters S n . 

The accuracy of the estimates for the phase shift parameters S n may be expressed 

20 as the variance of phase difference, as follows: 



\\ s {\<* nm -a m k\-ank) 2 w 2 (* > y)dxdy 



oU = 4^ „ , , : (26) 



wherein 
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= 5 7f 77 '/ — - (27) 



^ _ jj 5 Km {x,y)d xdy 

lis™ 2 ( x >y) dxd y 



Care is needed^n^evaluatiiig^ (26) and (27), as the phase y/'(x,y) is 

periodic and false discontinuities can disrupt the calculations if a nk near - k or k . 

Other methods for estimating these parameters c?^ that specifically incorporate 
5 the cyclic nature may be used instead of Equation (27), for example: 

&nk = { \\ s eX P [' Km " <*mk | ^ (*> >0 ^} ( 28 ) 

Fig. 4 shows a preferred embodiment of a automatic phase-shift calibration 
method 100. A first step 110 receives a sequence of fringe patterns f n . Each of the 

fringe patterns f n is real valued. An example of one of the fringe patterns f n from the 
10 sequence is shown in Fig. 5. Step 120 removes the fringe pattern offset a(x,y) from 
each of the fringe patterns f n . This is done by forming inter-frame differences g nm as 
defined in Equation (4). The inter- frame differences g nm are also real valued. 

Step 130 forms analytic images g nm corresponding to the inter- frame differences 
g nm as defined in Equation (15). This is preferably achieved by placing the inter- frame 
15 difference function g nm in the imaginary part of a complex image, and placing the 
imaginary output of the vortex operator V { g nm } in the real part. Fig. 6A and 6B show the 
magnitude and phase components respectively of an analytic image g nm of the example 
fringe pattern/, shown in Fig. 5. To enable the values of the phase components of the 
analytic image g nm to be displayed, the values have been adapted to a range [0;255] with 

127 

20 0 representing -n and 255 representing n . 

128 

Step 140 extracts the phase difference oc nm —a mk (modulo In) from two 
dependent analytic images g nm and g mk . Fig. 7 A shows a phase difference a nm - a mk 
between dependent analytic images g nm and g mk of the example fringe pattern f n . Again, 
the values of the phase difference a nm -a mk have been adapted to a range [0;255] with 0 
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representing -K and 255 representing n 9 allowing the values to be displayed. The 

128 

dark areas 710 and 711 have values close to zero on the representation, which corresponds 
to values for close to -n in the phase difference a nm -oc mk , whereas the light areas 720 
and 721 have values close to 255 in the representation, which corresponds to values close 
5 to 7t in the phase difference a nm - a mk . Accordingly, a value of 128 in the representation 
corresponds to a value of 0 in the phase difference a nm -oc mk . The spread of values in the 
representation is shown in the histogram of Fig. 7B. The peaks 730 and 731 in the 
histogram corresponds to phase difference a nm - a mk values. 

Step 150 determines a modulus of the phase difference \cc nm -a mk \. Fig. 8 A 
10 shows the modulus \ a nm -a mk \ of the phase difference a tm -oc mk represented in Fig. 7 A. 
Again the values have been adapted to a range [0;255] with 0 representing -/r and 255 
127 

representing tv . As can be seen from Fig. 8A, all the values have a small variation. 

128 

This is confirmed in the histogram shown in Fig. 8B, wherein the spread of values in the 
representation is shown. A single peak 810 in the histogram corresponds to modulus of 
15 phase difference \ a nm -oc mk \ values. It is noted that the values have only a small spread 

around the peak value. 

Step 160 calculates the estimated polarity factor // e (x,y), using Equation (21). 
The result is shown in Fig. 9 wherein all regions have values -1 or 1, with the dark areas 
corresponding to values of -1 and the light areas corresponding to values of 1 . The phase 
20 differences y ma are calculated using Equation (24) and the relative phase-shift parameters 

$n = $n-\ +Yn(n-\) 310 estimated in a consistent manner in step 170. 

The relative phase-shift parameters 5 n are then substituted into the general PSA 
defined by Equation (3) and the fringe pattern offset a (x,y) , the modulation b(x,y) , and 
phase y/'{x,y) are estimated in step 180. The resultant modulation b(x,y), phase 
25 i//'(x,y) and offset a(x,y) of the fringe pattern f n (x y y) of the example are shown in 
Figs. 10A, 10B and 11 respectively. It is noted that Figs. 10A and 11 appear almost 
uniformly white, which corresponds to almost no variation in the values of the modulation 
b(x,y) and offset a(x,y) with respect to coordinates (x,y) for the example shown. 
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Step 190 estimates the likely errors in the earlier contingent analytic images g am 
and g mk by means of a direct comparison of the analytic image g nm and the PSA image. 
Step 192 determines whether these errors are higher than a preset threshold. If the errors 
are smaller than the threshold, the method 100 ends in step 196. If these errors are higher 
5 than the threshold, the errors are incorporated into the weighting w(x,y) in step 194 and 
the method 100 returns to step 170 in an iterative manner, where the new weighting 
w(x,y) is applied to the estimation of the mean phase-difference y mn , until the errors are 
sufficiently low. In a simplest embodiment, the weighting w(x,y) is set to 1 for co- 
ordinates (x,y) where the comparison is good and 0 where the comparison is poor. 

10 Fig. 12 shows an example where the comparison of step 190 yields a weighting 

w(x,y) which is 1.0 where the comparison is good and 0.0 where the comparison is poor. 

It is generally noted that in equations with integrals over all values of x and y, such 
as Equations 5, 18, 19, 24 and 26, the integrals may be replaced by discrete summations 
for discretely sampled images. 

15 A general-purpose computer system 100, upon which the preferred embodiment of 

the invention may be practised, is shown in Fig. 13. The method 100 can be implemented 
as software, such as an application program executing within the computer system 200. 
In particular, the steps of the method are effected by instructions in the software that are 
carried out by the computer system 200. The software may be stored in a computer 

20 readable medium, including the storage devices described below, for example. The 
software is loaded into the computer system 200 from the computer readable medium, and 
then executed by the computer system 200. A computer readable medium having such 
software or computer program recorded on it is a computer program product. The use of 
the computer program product in the computer preferably effects an advantageous 

25 apparatus for estimating the phase of a sequence of fringe patterns, in accordance with the 
embodiment of the invention. 

The computer system 200 comprises a computer module 201, input devices such 
as a keyboard 202 and mouse 203, and output devices including a printer 215 and a 
display device 214. The computer module 201 typically includes at least one processor 

30 unit 205, a memory unit 206, for example formed from semiconductor random access 
memory (RAM) and read only memory (ROM), input/output (I/O) interfaces including a 
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video interface 207, an I/O interface for the printer device 215 and an I/O interface 213 
for the keyboard 202 and mouse 203. A storage device 209 is provided and typically 
includes a hard disk drive 210 and a floppy disk drive 211. A CD-ROM drive (not 
illustrated) may be provided as a non- volatile source of data. The components 205 to 213 
5 of the computer module 201, typically communicate via an interconnected bus 204 and in 
a manner which results in a conventional mode of operation of the computer system 200 
known to those in the relevant art. 

Typically, the application program of the preferred embodiment is resident on the 
hard disk drive 210, and read and controlled in its execution by the processor 205. 

10 Intermediate storage of the program may be accomplished using the semiconductor 
memory 206, possibly in concert with the hard disk drive 210. In some instances, the 
application program may be supplied to the user encoded on a CD-ROM or floppy disk 
and read via a CD-ROM drive (not illustrated) or floppy disk drive 211, or alternatively 
may be read by the user from the network (not illustrated) via the modem device (not 

15 illustrated). Still further, the software can also be loaded into the computer system 200 
from other computer readable medium including magnetic tape, a ROM or integrated 
circuit, a magneto-optical disk, a radio or infra-red transmission channel between the 
computer module 101 and another device, a computer readable card such as a PCMCIA 
card, and the Internet and Intranets including e-mail transmissions and information 

20 recorded on websites and the like. The foregoing is merely exemplary of relevant 
computer readable mediums. Other computer readable mediums may be practiced 
without departing from the scope and spirit of the invention. 

The method 100 of estimating the phase of a sequence of fringe patterns may 
alternatively be implemented in dedicated hardware such as one or more integrated 

25 circuits performing the functions or sub functions of estimating the phase of a sequence of 
fringe patterns. Such dedicated hardware may include graphic processors, digital signal 
processors, or one or more microprocessors and associated memories. 

Industrial Applicability 

30 The embodiment of the invention are applicable to interferometric testing of 

optical components and optical systems. The embodiment may also be used in non- 
destructive testing to study the deformation of objects under stress, and the study of 



504654.doc 



- 17- 



vibrating objects to identify vibration modes. Often, optical techniques are also used to 
visualise fluid flows, and variations in density are displayed in the form of fringe patterns. 

The foregoing describes only one embodiment of the present invention, and 
modifications and/or changes can be made thereto without departing from the scope and 
spirit of the invention, the embodiment being illustrative and not restrictive. 

In the context of this specification, the word "comprising" means "including 
principally but not necessarily solely" or "having" or "including" and not "consisting only 
of. Variations of the word comprising, such as "comprise" and "comprises" have 
corresponding meanings. 
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The claims defining the invention are as follows: 

1. A method of estimating relative phase shifts between fringe pattern images in a 

sequence of phase-related fringe patterns, said method comprising the steps of: 
5 a) removing offsets from each of said fringe pattern images to obtain pure 

AMFM patterns; 

b) determining contingent analytic images corresponding to each of said 
AMFM patterns; 

c) determining phase differences from dependent pairs of said contingent 
10 analytic images; and 

d) estimating phase shifts between pairs of said contingent analytic images, 
wherein said phase shifts between pairs of said contingent analytic images are said 
relative phase shifts between fringe pattern images. 



15 2. A method as claimed in claim 1 wherein said removing step is performed by 

calculating interframe difference images between pairs of said fringe pattern images. 

3. A method as claimed in claim 1 or 2 wherein said contingent analytic images 
comprise complex images with said AMFM patterns as imaginary parts and imaginary 

20 parts of a vortex operator applied to said AMFM patterns in said real parts. 

4. A method as claimed in any one of claims 1 to 3 further comprising the steps of: 

e) estimating a spatial phase of said fringe pattern images; 

f) estimating an amplitude modulation and said offset of said fringe pattern 

25 images; 

g) determining an accuracy measure of said contingent analytic images from 
said estimated amplitude modulation, offset and spatial phase; 

h) determining a weighting function wherein areas of low accuracy measure 
are given a low weighting; and 

30 i) iteratively repeating steps d) to h) by applying said weighting function in step 

d), until said accuracy measure is above a predetermined quantity. 
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5. A method as claimed in claim 4 wherein said accuracy measure determination 
step is performed by comparing one of said analytic images with a reconstructed fringe 
pattern, wherein said reconstructed fringe pattern is reconstructed using said estimated 
amplitude modulation, offset and spatial phase. • ,- r - - 

5 

6. A method of estimating a spatial phase of fringe pattern images in a sequence of 
fringe patterns, said method comprising the steps of: 

a) removing offsets from each of said fringe pattern images to obtain pure 
AMFM patterns; 

10 b) determining contingent analytic images corresponding to each of said 

AMFM patterns; 

c) determining phase differences from dependent pairs of said contingent 
analytic images; 

d) estimating phase shifts between pairs of said contingent analytic images; 

15 and 

e) estimating a spatial phase of said fringe pattern images. 

7. A method as claimed in claim 6 wherein said removing step is performed by 
calculating interframe difference images between pairs of said fringe pattern images. 

20 

8. A method as claimed in claim 6 or 7 wherein said contingent analytic images 
comprise complex images with said AMFM patterns as imaginary parts and imaginary 
parts of a vortex operator applied to said AMFM patterns in said real parts. 

25 9. A method as claimed in any one of claims 6 to 8 further comprising the steps of: 

f) estimating an amplitude modulation and said offset of said fringe pattern 

images; 

g) determining an accuracy measure of said contingent analytic images from 
said estimated amplitude modulation, offset and spatial phase; 

30 h) determining a weighting function wherein areas of low accuracy measure 

are given a low weighting; and 
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i) iteratively repeating steps d) to h) by applying said weighting function in step 
d), until said accuracy measure is above a predetermined quantity. 

10. A method as claimed in claim 9 wherein said accuracy measure determination 
5 step is performed by comparing one of said analytic images with a reconstructed fringe 

pattern, wherein said reconstructed fringe pattern is reconstructed using said estimated 
amplitude modulation, offset and spatial phase. 

11. A method as claimed in any one of claims 1 to 10 wherein an intensity of said 
10 fringe pattern images are in the form: 

L (*> y)=a(x,y) + b (x, y) cos \y (x, y) + S n ] 
wherein a(x,y) is said offset, b(x,y) is said amplitude modulation y/(x,y) is said 
spatial phase and 8 n is a phase shift parameter. 

15 12. Apparatus for estimating relative phase shifts between fringe pattern images in a 
sequence of phase-related fringe patterns, said apparatus comprising: 

means for removing offsets from each of said fringe pattern images to obtain 
pure AMFM patterns; 

means for determining contingent analytic images corresponding to each of said 
20 AMFM patterns; 

means for determining phase differences from dependent pairs of said contingent 
analytic images; and 4 

means for estimating phase shifts between pairs of said contingent analytic 
images, wherein said phase shifts between pairs of said contingent analytic images are 
25 said relative phase shifts between fringe pattern images. 

13. Apparatus as claimed in claim 12 wherein said means for removing offsets 
comprises means for calculating interframe difference images between pairs of said fringe 
pattern images. 

30 
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14. Apparatus as claimed in claim 12 or 13 wherein said contingent analytic images 
comprise complex images with said AMFM patterns as imaginary parts and imaginary 
parts of a vortex operator applied to said AMFM patterns in said real parts. 

15. Apparatus as claimed in any one of claims 12 to 14 further comprising: 
means for estimating a spatial phase of said fringe pattern images; 

means for estimating an amplitude modulation and an offset of said fringe 
pattern images; 

means for determining an accuracy measure of said contingent analytic images 
from said estimated amplitude modulation, offset and spatial phase; 

means for determining a weighting function wherein areas of low accuracy 
measure are given a low weighting; and 

iterative means for providing said weighting function to said means for 
estimating phase shifts, until said accuracy measure is above a predetermined quantity. 

16. Apparatus as claimed in claim 15 wherein said means for determining an 
accuracy measure comprises a comparator for comparing one of said analytic images with 
a reconstructed fringe pattern, wherein said reconstructed fringe pattern is reconstructed 
using said estimated amplitude modulation, offset and spatial phase. 

17. Apparatus for estimating a spatial phase of fringe pattern images in a sequence 
of fringe patterns, said apparatus comprising: 

means for removing offsets from each of said fringe pattern images to obtain 
pure AMFM patterns; 

means for determining contingent analytic images corresponding to each of said 
AMFM patterns; 

means for determining phase differences from dependent pairs of said contingent 
analytic images; 

means for estimating phase shifts between pairs of said contingent analytic 
images; and 

means for estimating a spatial phase of said fringe pattern images. 



504654.doc 



-22- 



18. Apparatus as claimed in claim 17 wherein said means for removing offsets 
comprises means for calculating interframe difference images between pairs of said fringe 
pattern images. 

5 19. Apparatus as claimed in claim 17 or 18 wherein said contingent analytic images 
comprise complex images with said AMFM patterns as imaginary parts and imaginary 
parts of a vortex operator applied to said AMFM patterns in said real parts. 

20. Apparatus as claimed in any one of claims 17 to 19 further comprising: 

10 means for estimating an amplitude modulation and an offset of said fringe 

pattern images; 

means for determining an accuracy measure of said contingent analytic images 
from said estimated amplitude modulation, offset and spatial phase; 

means for determining a weighting function wherein areas of low accuracy 
1 5 measure are given a low weighting; and 

iterative means for providing said weighting function to said means for 
estimating phase shifts, until said accuracy measure is above a predetermined quantity. 

21. Apparatus as claimed in claim 20 wherein said means for determining an 
20 accuracy measure comprises a comparator for comparing one of said analytic images with 

a reconstructed fringe pattern, wherein said reconstructed fringe pattern is reconstructed 
using said estimated amplitude modulation, offset and spatial phase. 

22. Apparatus as claimed in any one of claims 12 to 21 wherein said fringe pattern 
25 images are in the form: 

/„ (x, y) = a (x, y ) + b (x, y) cos [y/ (x, y) + S n ] 
wherein a(x 9 y) is said offset, b(x,y) is said amplitude modulation y/(x, y) is said 
spatial phase and 8 n is a phase shift parameter. 

30 23. A method of estimating relative phase shifts between fringe pattern images in a 
sequence of phase-related fringe patterns, said method being substantially as described 
herein with reference to the accompanying drawings. 
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24. Apparatus for estimating relative phase shifts between fringe pattern images in a 
sequence of phase-related fringe patterns, said apparatus being substantially as described 
herein with reference to the accompanying drawings. 

Dated 2 August, 2000 
Canon Kabushiki Kaisha 

Patent Attorneys for the Applicant/Nominated Person 
SPRUSON & FERGUSON 
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